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Abstract 

We apply both a scalar field theory and a recently developed transfer-matrix 
method to study the stationary properties of metastability in a two-state 
model with weak, long-range interactions: the Nxoo quasi-one-dimensional 
Ising model. Using the field theory, we find the analytic continuation / of the 
free energy across the first-order transition, assuming that the system escapes 
the metastable state by nucleation of noninteracting droplets. We find that 
corrections to the field-dependence are substantial, and by solving the Euler- 
Lagrange equation for the model numerically, we have verified the form of the 
free-energy cost of nucleation, including the first correction. In the transfer- 
matrix method we associate with subdominant eigenvectors of the transfer 
matrix a complex-valued "constrained" free-energy density computed di- 
rectly from the matrix. For the eigenvector with an associated magnetization 
most strongly opposed to the applied magnetic field, fa exhibits finite-range 
scaling behavior in agreement with / over a wide range of temperatures and 
fields, extending nearly to the classical spinodal. Some implications of these 
results for numerical studies of metastability are discussed. 

PACS Number(s): 64.60.My 64.60.Qb, 02.70.Rw, 03.50.Kk 
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I. INTRODUCTION 



Determining the stationary properties of metastable states from the standpoint of sta- 
tistical mechanics has been the topic of many studies over the last three decades. (For 
reviews, see e.g. Refs. IHJ^.) Some early treatments of this problem have shown that tra- 
ditional methods of equilibrium statistical mechanics in a more generalized form might be 
applicable to metastable states. In a study of the analytic properties of the free energy at 
the condensation point, Langer 0] conjectured that the imaginary part of the free energy 
F analytically continued from the equilibrium phase across the first-order phase transition 
may be associated with the decay rate of the metastable phase. A dynamical investigation 
of thermally activated nucleation showed that for a wide class of models, the decay rate 
r may be written as 

r = ^|ImF|, (1) 

TT 

where /3 is the inverse temperature l/k-QT, and k is a kinetic prefactor that contains all 
dependence on the dynamics. 

Subsequently there have been a number of theoretical and numerical studies of metastable 
decay. Binder and collaborators developed a scaling theory based on a proposed defi- 
nition of metastable states in terms of a nonequilibrium relaxation function and tested this 
theory by Monte Carlo simulation on the 2D Ising model. Schulman and coworkers studied 
by various methods metastability in the 2D Ising model in the ID Kac model with alge- 
braically decaying interactions |jlO| , in the Curie- Weiss model , and in a dropletlike "urn" 



model |jT2[ . In addition, Biittiker and Landauer ITSlJT^ studied nucleation in the overdamped 
ID sine-Gordon chain, and Klein and linger [15,16| used a 0^ field theory to study systems 



with long-range interactions in arbitrary dimensions near the classical spinodal. Other stud- 
ies of nucleation in short-range systems also exist in the literature, including field-theoretical 
studies [^^-|19[, series expansions [pO|-p3|, and exact diagonalization studies |^. Each of the 



above studies supported the validity of hanger's treatment. Recently Gaveau and Schulman 
23| have determined a rigorous upper bound for F for a larger class of models than that 



considered by Langer, and have used it not only to explain why Eq. (|l]) is usually valid, but 
also to provide an example in which it can give misleading results. 

A study of the analytic properties of transfer-matrix eigenvalues led to a conjecture 
p6| that there is a correspondence between the analytically continued free-energy density 



/ and the analytic continuation of the dominant eigenvalue of the transfer matrix. This 



conjecture has been supported by subsequent numerical studies [Pi p7| , p8| of the 2D Ising 
model. In a similar study of an Ising model with Nxoo cylindrical geometry in which the 
interaction range is linear in the cross section, evidence from the transfer-matrix eigenvalue 
spectrum was found for the emergence of a classical spinodal in the limit of weak, long- 
range interactions. In the present work we apply a recently developed "constrained-transfer- 
matrix" (CTM) method ||30H33| to study the properties of the metastable phase of this quasi- 
one-dimensional Ising (QIDI) model. Since the transfer matrix for the Nxoo QIDI model 
has a rank of [^, instead of the typical value of 2^ for short-range Ising systems, 

it is relatively easy to study large systems by transfer-matrix techniques. Using the CTM 
method, we obtain a "constrained" free-energy density associated with the metastable phase. 
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Our main focus here will be to relate the finite-size (or equivalently for this model, finite- 
range) scaling behavior of the constrained free energy to the scaling behavior of the analytic 
continuation of the free energy for an equivalent one- dimensional long-range field-theoretical 
model. This work represents an extension to metastable states of our recent study ||3^ of 
critical equilibrium finite-range scaling in the QIDI model. 

In studies of condensation it is often useful to consider the behavior of long-range models. 
One reason for this is that long-range models share many of the equilibrium and metastable 
properties of mean-field models, which are often exactly soluble. In addition, as with equi- 
librium states near the critical point, the metastable states of these models near the classical 
spinodal may be studied using a simple field-theoretical Hamiltonian [|I5],|T6|. Further, the 
behavior of many physical systems, such as coherent metal- hydrogen systems ||35|, superfiuid 
^He [^, long-chain polymer mixtures p7| , |38[] , materials exhibiting elastic phase transitions 
39| , superconductors [^, and ferroelectrics with long-range strain fields p] , ^ , indeed are 
well described by long-range models with suitably chosen order parameters. Many treat- 
ments of classical metastability for these models are reported in the literature, using many 
different approaches: the Fokker-Planck equation p!l|J^J^ , the properties of the nonequi- 



librium relaxation function 



45 



simulation 



and Monte Carlo 

Tn,ira,ira 



renormalization group analysis 
in addition to analytic continuation of the free energy 
The approach of the present work is to calculate numerically from the transfer matrix a 
constrained free-energy density 



,31 



associated with the classical metastable state, and 
to determine its scaling behavior in weak fields away from the critical point, and near the 
classical spinodal. We demonstrate that this quantity is closely related to the analytically 
continued free-energy density. 

The remainder of this paper is outlined as follows. In Sec. |I| we derive the Euler- 
Lagrange equation for the order parameter in the QIDI model, and we identify the solutions 
corresponding to the equilibrium and metastable phases, as well as the most likely fiuctuation 
through which a QIDI system in the metastable phase can decay into the equilibrium phase. 
In Sec. |ITl| we obtain a mapping of the QIDI model near the classical spinodal to a one- 
dimensional (fr" field theory. Using both the solution to the Euler-Lagrange equation and the 
field theory, we perform the analytic continuation of the free energy for the model across 
the first-order phase transition, obtaining an exact expression for its imaginary part near 
the classical spinodal. In Sec. ^ we give a brief overview of the CTM method, which is 
based on the notion of "constrained" joint probability densities [pO| , ^ , and we show how 
the symmetries of the QIDI model |]2S| are used to simplify the calculation. In Sec. |V] we 
verify by numerical integration of the Euler-Lagrange equation the temperature and field 
dependence of the the free-energy cost of nucleation, including the first field correction. 
We also present our numerical transfer-matrix results and compare them with finite-range- 
scaling predictions for / based on the field-theoretical results of Sec. |T|. In Sec. ^ we 



summarize our conclusions, and we discuss the implications of the results to the applicability 
of the CTM method to studies of metastability. 



II. THE QIDI MODEL 



The quasi-one-dimensional Ising, or QIDI, model ||2^ is a one-dimensional chain of L sub- 
systems, or layers, each of which contains Ising spins Si „=±1, where the index i=l, . . . ,L 
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runs over the layers, and the index n=l, . . . ,N runs over the spins in a given layer. Each 
spin interacts ferro magnetically with each of the 2N spins in the adjacent layers with a 
coupling J/N>0 and with the other N—1 spins of its own layer with a coupling J'/N. Each 
spin also interacts with an external magnetic field H. The Hamiltonian is explicitly 

L / J N N ji N-l N N \ 

1=1 \-^^ n=l n'=l -'^ n=l n'=n+l n=l / 

where periodic boundary conditions are imposed, that is (Wi,n) SL+i,n=Si,n- The sums over 
n and n' in Eq. (|^) may be performed directly, so this Hamiltonian may be expressed in 
terms of a discrete field of layer magnetization densities 0j=A^~^ Y.n=i ^i,n as 

L 

'H[(t>] = -NY. [J<^>^<^>^+l + \J\<PI - N~^) + HcP^ . (3) 
i=l 

The partition function may then be written as 

Z = Y.9W-'''^''K (4) 

where 

*Hn(.^,f,J (5) 

denotes the multiplicity of spin configurations giving rise to 0. In the limit of large A^, the 
probability density of is sharply peaked around its maximum, so if we define a free-energy- 
density functional as 



^[0] = {NL)'^ (7^[0] - ln(^[0])) , (6) 
then the free-energy density is 

/qidi = -{l3NLr^ \nZ^ min^[0] . (7) 

Applying Stirling's approximation to Eq. (H), combining the result with Eqs. (H) and (H), 
and removing constant terms, we have for the QIDI model 



•^M = ^ E [\J{<P^+l - 00' + /CW(0.)] , (8) 



i=l 

where 

M,) = - i/*. + r ' (i±* In i±* + 1^ 1„ i^) (9) 

is the free-energy density of a Curie- Weiss ferromagnet with a critical temperature 

= 2 J + J'. (10) 
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Corrections to few due to finite will be briefiy considered in Sec. 

In Ref. [0] we showed, both analytically and numerically, that the critical equilibrium 
iV-scaling relations for the QIDI model are those of an Nxoo cylindrical Ising model above 
its upper critical dimension, which can be obtained from a standard one- dimensional 0"^ field 
theory by a single rescaling of the correlation length. Here we extend our study of the scaling 
behavior to metastable states, where P>Pc and Hj^O, by studying the scaling behavior of the 
analytic continuation of the free energy across the first-order phase transition. To do this, 
we must first calculate the free-energy density for the equilibrium and metastable phases of 
the model, as well as the height of the free-energy barrier between them. For large N, we 
can treat each 0j as a continuous variable. In this case, extremizing J-'[4>] with respect to 
in the L-cube [—1, 1]^, we obtain a system of L coupled nonlinear equations: 



+ 



+ 



df 



cwi 







(11) 



We can gain some understanding of the solutions to Eq. (|n|) by expressing the equation in 
differential form. We define a function (j){r) to be continuous in a dimensionless longitudinal 
coordinate r, and we force (j){r)=(j)i at integer values r=i. By Taylor's theorem, if z < r < 
i + 1, then — (piY = (V0(r))^ + O(V0V^0). Therefore, if 0(r) does not vary too 

rapidly over the length of the chain, then the free-energy density is well approximated by 
the Ginzburg-Landau- Wilson form: 



1 
L 



dr 



i J(V0)2 + /cw(0(r)) 



which represents an effective interaction range 



R 



N 



(12) 



(13) 



as measured in sites. Note, however, that Rn is finite as measured in layers. In the critical 



region, this mapping to the continuum field theory agrees with the formulation in Ref. [34 



The resulting Euler-Lagrange equation for stationary solutions is the continuous analogue 
to Eq. (inD: 



6(f){r) dr"^ 
which can be integrated once to give 







(14) 



(r) 



dr 



fcw{4>{r))-C 



(15) 



where C is a constant. 

Equation (|T^ may be interpreted as the equation of motion for a particle of mass J 
moving in a potential —fcw{(p{^))y where 0(r) represents the displacement of the particle 
at time r. Detailed descriptions of the solution from this point of view can be found in 
Refs. In determining the solutions, we assume that 0(r) is continuous. Since few is 

continuous on the interval [—1,1] and continuously differentiable on (—1,1), 0(r) is twice 
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continuously differentiable on (—1, 1). The boundedness of requires that if the LHS of Eq. 
(p!5|) is nonzero for all r, it must approach zero as |r|^oo, since would in that case be 
monotone, while the nonnegativity of the LHS of Eq. (|I5|) requires that (Vr) /cw(0('^))>C- 
The first condition on the solution is therefore that the range of is an interval /c[— 1, 1] 
for which inf (/cw(-^))=C- Further, we must have as a second condition that fcw{(p)=C 
at the endpoints of /. To see this, assume on the contrary that fcw{(p')>C for one of the 
endpoints 0'. For some r we would then have (f){r)=(f)' and \d(f) / dr\>0 by Eq. (PTSP, forcing 
0(r) to violate the first condition by taking values outside the interval /. 

By varying the constant C we find three types of solutions satisfying the above con- 
straints. These are illustrated in Fig. |I|. The first type (I) is a completely uniform solution, 
in which dcf)/ dr=d?(f)/ dr'^=Q for all r. By Eq. (|T^ , these solutions exist only if the constant 
value of 0(r) is an extremum of /cw- There can be as many as three solutions of this type, 
and they are marked on the free-energy sketches in Fig. |I| as diamonds. The second type of 
solution (II) is a nonuniform "droplet" or "interface" solution, in which one of the endpoints 
of / is the locus of a local minimum 0ins of /cw that is not a unique global minimum. If such 
a minimum exists, there is only one solution of this type. The vertical lines descending from 
the free-energy sketches in Fig. |l| mark the endpoints of the range of 0(r) for this solution, 
and sketches of the solution in real space are drawn below. If the type-II solution exists, 
then a band of solutions of a third type (III) also exists. These solutions are oscillatory in 
space and are marked on the free-energy sketches in Fig. |I] as hashed regions. They are not 
considered in the following discussion since the free-energy densities associated with these 
solutions are the highest of those for the three types. 

By considering the solutions of types I and II, we easily identify the equilibrium and 
metastable states as those with the two lowest values of JF. Both are uniform (type-I) 
solutions, 0(r)=0eq and 0(r)=0ms, respectively, and thus by Eq. (|TT|) have ^[0]=/cw('/'eq) 
and JF[0]=/cw(0ms), respectively. The metastable state exists at temperatures and fields 
for which /cw exhibits two minima. For (3>l3c, this is true if 0ms does not coincide with a 



point of inflection of few- The locus of the inflection is given by = iyl — (3c/ f3, so the 
metastable state exists if 0^^ > 1 — /3c//3. Combining this constraint with Eq. (p!T| ) and with 



the uniformity of 0, we obtain the classical (Curie- Weiss) spinodal field ||rT],g4|, which 
defines the limit of metastability: 

= -I3;'<j), + (3-' tanh-i 0, , (16) 

where 



0s = -(sgni/)^l-/5,//3. (17) 

The solution with the next-lowest value of is the type-II solution. This solution represents 
the fluctuation of lowest free energy through which the metastable phase can decay into the 
equilibrium phase (for Hy^O), or through which a system can pass between two coexisting 
phases (for H=0). 

For 0<|if|<|ifs|, we can estimate the free-energy cost AF of this fluctuation by numer- 
ically solving the coupled nonlinear equations defined by Eq. ([TT|) for an L-layer system, 
allowing the magnetization densities to vary continuously. In Sec. [V| the numerical integra- 
tion procedure is outlined and its results for AF are compared with the analytic result from 
a field-theoretical Hamiltonian near the spinodal and with our transfer-matrix results. 
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III. ANALYTIC CONTINUATION OF F 



For the QIDI model with /3>/3c, the free energy F in the hmit N^oo is not everywhere 
analytic, but rather exhibits a discontinuous first derivative with respect to H at H=0. 
However, if we vary H continuously through H=0, an analytic continuation / of the free- 
energy density across H=0 does exist as an analogue of Eq. (^, where the partition function 
is constrained to configurations that do not allow the system to reach equilibrium. This 
continuation is the minimum of that coincides with /qidi at H=0, but increases as H 
departs from zero. When H reaches the spinodal field, the metastable minimum vanishes. 
At this point the continuation has a branch point and becomes complex as \H\ is increased 
further. 

We determine the leading behavior of the analytic continuation of F in the region of the 
spinodal {(j)=(f)s, H=Hs) as follows. Taking the continuum limit used in the previous section, 
we write 0(r) = (j)s+v{r) and H = H^ + X, and we expand Eq. for the free-energy-density 
functional near the spinodal to third order: 



J^[T,H,^]=J^[T,H,,4>,]+AJ^ 



(18) 



with 



1 

AJ^ = - / dr 
L Jo 



U{\/vf + X 



cw 



OH 



+ Xv 



cw 



d<pdH 



1 

-A0S + T / '^^ 

L Jo 



where we find 



a 



/3 



'1 - 



k2\2 



/5; 



/5 



(19) 



(20) 



using the explicit form of Eq. (|^) for /cw and Eq. (|T^. In Eq. (|T^) we have shown only 
the derivatives that are not identically zero. (Since /cw is stationary and has an inflection 
at (f)=4>s, the first and second 0-derivatives are among those not shown.) From the above 
expansion, the requirement that is stationary gives the Euler-Lagrange equation: 



- JV\ -X + av^ + 0{v^) = 



(21) 



Since the critical fluctuation consists of a single droplet, the free-energy cost of the 
fluctuation is not extensive in L. We take the position of the droplet core to be r=0 in 
the following discussion. Since this fluctuation is local, the solution to Eq. (|2lD must be 
asymptotically uniform in the limit L^oo, that is f (r)— ^constant as |r|^oo. From Eq. ( pi]) 
we thus have asymptotically 



v^vo = {X/ay/^ + 0{X) 



(22) 



with the sign restriction sgnfo=sgna if |//|<|ifs|- Changing variables to u{r) = v{r) — vq, 
the Euler-Lagrange equation becomes 
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JV^M + 2avou + av? + 0{v?) = . 



(23) 



In the limit |r|— i>oo we may neglect terms of O(-u^). We thus obtain 

n(r) ~ exp(-r/^r) , 

where 

6 = 



|aA|-^/^ + 0(Ai/^) 



(24) 



(25) 



is the relaxation length of the fluctuation. The power-law dependence on A of the order 
parameter vq and the length scale is indicative of a critical point at the spinodal, as was 
pointed out previously in Refs. []T3,TB,|3|. We can use the change of variables above to 
separate the L-extensive part of the integral in Eq. (p!9|): 



-A0S + - dr I J(Vm)2 + a{-lvl + v^u^ + \u^) + 0{v^) 
L/ Jo 



1 
L 



dr 



(26) 



where we have integrated (Vm)^ by parts and used Eqs. (p2|) and (^). 

If |if|>|ifs|, then vq is purely imaginary. For the uniform solutions of the Euler-Lagrange 
equation, f =±i|fo|, AjF has an imaginary part given, to this order of approximation, by the 
second term in Eq. (pGl) . Combining this with Eq. ( ^21) and noting that Vq is real, we have 

(27) 



Imf = Im(A^) = ±||ari/2|^|3/2 + 0{X'/^) 



This expression is exactly the result for a Curie- Weiss mean-field model [11,44 



On the other hand, if |if|<|ifs|, then AjF is purely real, and the first three terms of Eq. 
(^) give the contribution of the uniform metastable background, whereas the integral term 
gives the contribution of the critical droplet. Taking the limit L—>-oo, we insert the explicit 



solution to Eq. (|2l|) |14HT6 



fi(r) = t>o 1 — 3sech^(r/2^i.) 

into the integrand in Eq. (^) to find the free-energy cost AF of forming the droplet: 

AF = fN^^\a\~'/'\X\'/^ [l + 0{X'/^ 

/2J|a|-=^/^|A|5/^ + 0(A^/^) . 



(28) 
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(29) 



The first expression for AF in Eq. (^) illustrates the difference in the behavior of the 
analytically continued free energy between long-range models, for which the length scale for 
the critical droplet is a characteristic length ^r'^L ||T5|JT6[] , and mean- field models, for which 
the only length scale for fluctuations is the length L of the entire system. The volume- 



extensivity and field-dependence of the mean-field result for AF [p!l|^ can be recovered by 
replacing ^j. with L in Eq. (PU]). 

Note that if |/7|<|ifs|, the analytically continued free energy F is real- valued only in 
the limit of infinite interaction range. However, we are interested in the scaling of the free 
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energy for systems with a long, but finite interaction range. For such systems F is complex, 
but its imaginary part approaches zero rapidly as Rn—>-oo. We can determine the scaling 
properties of the free energy by expanding the partition function for the system under the 
constraint that all fluctuations remain subcritical. Following Refs. P,^, we write 



Z = Zo + Zi + Z' 



(30) 



where Zo is the contribution from the region of configuration space around the local minimum 
of corresponding to the metastable background, Zi is the contribution from the region 
around the saddle point corresponding to the critical droplet, and Z' represents higher- 
order terms. The free energy must be extensive in L, so we must also consider contributions 
of multiple critical droplets. If we assume that Z' has contributions only from multiple 
identical, noninteracting droplets, then by expanding Z as a series in Zi/Zq we have ||4^ 



Z = Zoexp(Zi/Zo 



(31) 



For the stationary points of JF identified above as the uniform metastable background Vq 
and the critical droplet Vi, we write v{r) = f„(r) + z/(r) and expand to second order: 



(32) 



where 



L 



dr 



1 2 9'^ few 



- / dr UiVuf + (avn + O(f^) 
L Jo ^ ^ 



(33) 



In Eq. (|33|) we have expanded the second term of the integrand to first order in u„ since 
the values of Vn{r) are small, and we have again used the fact that few has an inflection at 
4>=4>s- We diagonalize the quadratic form by a principal-axis transformation i'{r)=J2j CLj^j 
to the orthonormal set of eigenmodes of the Schrodinger equation 



(34) 



From Eqs. (^2])-(0) the resulting contribution Z„ to the partition function can thus be 
written as a product of decoupled Gaussian integrals: 



dvexp{-pNLJ^[v]) 



exp{—/3NLT[vn\)Y\_ J daj exp{—PNa'^ujj 



(35) 



Note that although the results that appear later in this section are restricted to the region 
of the spinodal, Eq. (^) holds everywhere in the region of metastability. It can easily be 
seen from this and from Eq. (31) that in this region a Boltzmann weight equal to e~^^^ 
appears in the analytically continued free energy. 

In order to determine the behavior of the analytically continued free energy in the region 
of the spinodal, where AF is small, we must further determine the eigenvalues of Eq. (p4|). 
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For n=0 Eq. (|3^) is a free-particle equation with a potential floor at avo, whereas for n=l 
it describes a particle in a potential well described by avi{r). There are three localized 
eigenmodes for n=l |T^. They are 



uq oc sech^(r/2^r) 

z^i oc sech^(r/2^r)tanh(r/2,^r) 

U2 oc 4sech(r/2^r) - 5sech^(r/2^r 



(36) 



for which 



UJo 
UJ2 



+ O(A) 



-l\aX\'/^ + 0{X) 
\aX\^/^ + 0{X) 



(37) 



are the respective eigenvalues [4£]. The first eigenmode uq is the only unstable mode of 



fiuctuation and corresponds to an increasing or decreasing magnetization at the droplet core. 
The second eigenmode i^i is a translational mode, and the third eigenmode 1^2 corresponds 
to a widening or narrowing of the droplet. All other eigenmodes are unbounded and form a 
continuous spectrum with Wj>iJ/^^- 

We shall consider the contributions of each localized eigenmode individually. For each 
ujj>0, the Gaussian integral of Eq. ( pSf ) is a well-defined function, G{ujj) = {n/PNujjY^'^, but 
if u)j<0, as is the case for j=0 and j=l, the Gaussian integral is divergent. We handle each 
divergence in a different way. For j=0, we analytically continue G{uj) from the half-plane 
defined by Rea;>0 to the entire complex plane. The result is purely imaginary and has a 
sign ambiguity, since the continuation of G has a branch cut on the negative real axis. For 
j=l, we use the fact that ui is the normalized derivative of Vi (Eq. (^)) with respect to r. 
We can thus make a transformation of the integral to one over dr, the Jacobian of which is 
simply ||Vfi||. The resulting contributions to Zi from the localized modes are 



± 



TC 



2^ f3N\uJo\ 



(j = 0) 



lAF 



L (J = l) 



PNU2 



(j = 2) . 



(38) 



The factor of | has been introduced in the j=0 result |^ to account for the fact that the 
integral inappropriately counts a divergence as one moves toward the metastable minimum 
48| . The form of the j 
argument [p!7|,p!S 

spectra of eigenmodes. Using a previous calculation of the ratio of the relevant eigenvalue 
products in App. B of Ref. |T^, we obtain 



1 result ||50| is quite general and follows from a simple virial 
We are then left with the contributions to Zi/Zo from the continuous 



30J\uJo\uj2 ( (3N^ 



3/2 



n 



00 



6 



TT 



(39) 



10 



where tu^ denotes an eigenvalue of the free-particle Schrodinger equation. Combining this 

, and (^8|), we have 




2np N^, 
±12{V2Jnl3N)-^^'^\ 



a 



-1/81 M 7/8 



1 + 0(A 



l/2^ 



-I3AF 



(40) 



Assuming that the dynamics of the system is governed by a Fokker-Planck equation, we 
have from Ref. that the kinetic prefactor K=f5'y\(jJo\, where 7 is the fundamental rate of 
fluctuation. With Eqs. (|1|), (0), and (|iO|), this gives the nucleation rate density 



NL 



157(v^iV)-^/^ (f ) 



3/2 



\a 



3/81 X 111/8 



|A|^^/« 1 + 0(A 



1/2N 



-/3AF 



(41) 



As was pointed out earlier, for large N, unless H is extremely close to Hs, the free-energy 
cost AF of surmounting the nucleation barrier is large, so the exponential factor sets the 
scale for the metastable lifetime. However, for small N, or for H^Hg, the lifetime is more 
strongly dependent on the dynamics and on the detailed structure of the saddle point. In 
Sec. we will show how the crossover between these two regimes depends on by deriving 
the finite-range-scaling properties of Im/. In the same section we will also directly compare 
numerical results from the transfer-matrix method outlined in the next section with the 
exponential weight in Eq. (^OD, where the gap AF is determined numerically from Eq. (|11]). 



IV. THE CONSTRAINED-TRANSFER-MATRIX METHOD 

Recently one of us |l30| , |3T| introduced a transfer-matrix method, based on the concept 
of "constrained" joint probability densities, to obtain analogues of the free-energy density 
for constrained states. Preliminary applications |l30|- p^] to the QIDI model have shown 



qualitative agreement between the behavior of the free-energy- density analogue associated 
with the metastable eigenvalue branch of the transfer matrix and the analytically continued 
free-energy density. In this study we use this constrained-transfer-matrix (CTM) method 
with finite-range scaling |Q to obtain more quantitative results for the scaling of the imag- 
inary part of the metastable "constrained" free-energy density. In this section we outline 
the method, and we describe the use of symmetry reductions in its application to the QIDI 
model. 



The standard transfer-matrix method is usually applied as follows |£T|. Starting with 
a Hamiltonian for an NxL g-state system that is invariant under translation in the L- 
direction, such as Eq. (0), we write it as a sum of layer Hamiltonians: 'H=Y,f^i'H{xi, Xi+i). 
The layer Hamiltonians Ti. are chosen to depend only on the configurations of a pair of 
adjacent layers and are chosen so that the form of Ti, is independent of the layer index. We 
then define the q'^ xq^ transfer matrix T as an operator in the dual space of configurations 
\xi) = \si^i) ■ ■ ■ \si^N) of each of the two adjacent layers: 

T = ^|x)e-^^("'"')(x'| . (42) 
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The partition function for the L-layer system with periodic boundary conditions is then 
Z=Ty{T^), which in the hmit L— >-oo gives the free-energy density in terms of the largest 
eigenvalue Aq of T as /=— (A^/?)~^ In Aq. (Here we index eigenvectors a=0, 1,2, . . ., in 
order of decreasing magnitude of their respective eigenvalues: Ao>|Ai|>|A2|>- ■ ■•) Since T 
is positive, Aq is positive and nondegenerate, and the eigenvector |0) can be chosen to have 
all positive elements by the Perron- Frobenius theorem The equilibrium state of the 

system is characterized by the joint and marginal probability densities 

Poixi,Xi+k) = (0|xi)(xi|(Ao^T)l^l|a;i+fc)(a;i+fe|0) 

Poix,) = (0|a;,)(x,|0) . (43) 

In the following discussion we restrict T to be symmetric. (For the QIDI model T is 
made symmetric by symmetrizing the layer Hamiltonian 7i.) For each eigenvalue A^ of T 
we define a constrained transfer matrix (CTM) T^, to commute with T, so that it can be 
expanded in the eigenvectors \a') of T: 

T^ = Y.\a')fi^{a'){a'\ . (44) 

a' 

The "reweighted" eigenvalues are chosen to produce "constrained" probability densities, 

Pa{xi,Xi+k) = (a|xi)(xi|(A^^Ta)l^l|xi+A:)(a;i+fc|a) 

Pa{xi) = {a\xi){xi\a) , (45) 

in analogy with the equilibrium (a=0) case. It was pointed out in Refs. P J27| , |28[] that 
the constrained marginal probability densities Pa{x), as defined above, can be interpreted 
as actual probability densities of single-layer configurations in a constrained state. For 
example, the expectation value (M)^ of the layer magnetization density for a constrained 
state I a) is given by 

{M)^ = Y,P^{x)M{x) . (46) 

X 

In order to ensure that the entire system is characterized by Pa{x), the matrix Tq, 
must be chosen so that the constrained joint probability densities Pa{x,x') satisfy the fol- 
lowing regularity conditions: (^) that Pa{x) can be obtained by summing over the config- 
urations of one layer, Pa{x)=J2x' Pa{x,x'); (ii) that Pa{xi,Xi-^k) is well-defined for k=0, 
Pa{xi,x'j)=6x-^x'-Pa(yXi); and (m) that Pa{xi,Xi-^k) reflects stochastic independence in the 
limit s>oo, \im\k\^oo Pa{xi, Xi^k)=Pa{xi)Pa{xi+k) ■ For a matrix Tq, chosen to commute 
with T, these requirements are satisfied if has the same rank as T, and if its dominant 
eigenvalue is Aq,. The reweighting scheme used in this work is 

Uc^') = I ^'/^"' S^' < M (47) 

This choice ensures that any fluctuation described by Pa'{x), with a'y^a, decays away to 
the |a) eigenstate with a length scale | ln(AQ,/AQ,/)|~^. Thus a fluctuation with a'<a, which 
would otherwise grow, is effectively suppressed. We note, however, that the choice of T„ 
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given by Eqs. ( PD and (^T)), which satisfies the constraints, is not unique. It is easy to see 
that we recover the original transfer matrix T when we set a=0. 
The "constrained" free-energy density fa is defined by 

f. = {n)a-f3-'Sa, (48) 

where 

{'H)a = ^ J2 Pa{Xu Xi+i)n{Xi, Xi+i) (49) 

is the expectation value of H with respect to Pa{x,x'), and a generalized entropy Sa is 
defined using Pa{x,x') in analogy with the source entropy of a stationary, ergodic Markov 
chain (See e.g. Ref. |]53|): 

Sa = -^ Pa{Xi,Xi+i)Ln{Xi\\a^Ta\Xi+i) . (50) 

The generalized free-energy density may be written in the form 

ln|AQ,| 1 (xi|TQ,|xj+i) 

The first term is analogous to the equilibrium case, and the second term is complex-valued 
in general. This is because is not a positive matrix in general. To see this, note that 
for |a) (07^0) to be orthogonal to |0), the elements of |a) must be of mixed sign. Since 
is the largest eigenvalue of Tq,, the largest contribution to Tq, is the projection |a)AQ,(a|, 
which must contain negative elements. Therefore, the argument to the principal value of the 
logarithm (Ln) may be negative. We define the domain of Ln^; to be |2;|>0, —7c<aig z<7r, 
thus choosing the branch cut along the negative real axis. It is easy to see that if a=0, the 
second term of Eq. (^TJ) vanishes, leaving the equilibrium free-energy density. 

The above formalism is applicable to any system for which a symmetric transfer matrix 
can be written. A CTM study of the two-dimensional Ising ferromagnet with nearest- 
neighbor interactions is reported in Ref. |^3[, and a similar study of a three-state model 



with long-range interactions is reported in Ref. One reason we have chosen to study 
the QIDI model is that the Hamiltonian is invariant under any permutation of spins in a 



given layer, and thus the rank of the transfer matrix for this model is p9|], instead of 

the typical value of 2^ for short-range Ising systems. We use the low rank of the transfer 
matrix to our advantage by expressing the matrix in a reduced basis. The computer time 
and memory saved by this reduction allow us to study systems with very large cross sections. 
The reduction proceeds as follows. Let X be the 2^-dimensional vector space in the basis 
(i=l) • • • ) 2^) of layer configurations, and let G be the group of spin permutations 
of the layer, which are represented by unitary operators U^-.X^X (k^G). This group 
partitions into equivalence classes C„ (^=0, . . . , N), each of which may naturally 

be associated with a magnetization m„. This partition induces a decomposition of the 
configuration space as X=0X„, where each X„ is the subspace spanned by the basis vectors 
|xj)GC„. That T is G-invariant means that (VfcGG) TUjfc=T, which implies the weaker 
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property that T commutes with every Ufc, and thus may be simultaneously diagonalized 
with any Ufc by a unitary transformation ■ S^^. Thus if Dx and are diagonal matrices 
containing the eigenvalues of T and respectively, then DTD;j=SfcTUfcS^^=SfcTS^^=DT. 
Since the eigenvalues of Ufc are roots of unity, it follows that for each eigenvalue of U/j 
not equal to unity, the corresponding transfer-matrix eigenvalue is zero. Therefore, T has 
nonzero eigenvalues only in the subspace of X for which the eigenvalues of Ufc are unity for 
every /cgG. This is the A^+l-dimensional G-invariant subspace of X. 

We construct a basis in this subspace as follows. For each class C„, let |m„) denote the 
normalized G-invariant vector projecting only into X„: 

\mn) = Qn"^ E l^i) ' (52) 

where gn={^ is the number of basis vectors in C„. Using Eqs. ( ^2] ) and (|5^), we rewrite 
the transfer matrix in this reduced basis as 

T = E E E k)e-^^(^'^')(x'| 

n,n' x£Cn x'eC^i 

= E l^n)(^7n^?n')'/'e-^^(™"'"^"')(m„.| . (53) 



We can similarly decompose the sums over configurations present in Eqs. (^) and (|49|-51) 



into sums over the reduced basis vectors [0. The transfer-matrix results presented in the 
next section were obtained using this reduced representation. 



V. NUMERICAL RESULTS 

In this section we present the results of two numerical methods used to study the prop- 
erties of the metastable phase in the QIDI model. We numerically integrate the Euler- 
Lagrange equation, in the discrete form of Eq. (|Tl|), to estimate the free-energy cost of 



nucleation, which we use to verify the field-theoretical results of Sec. pT[ We also apply 
the CTM formalism outlined in Sec. ^ to the model. We perform finite-range scaling on 
the constrained free-energy density fa of Eq. (|5lD , and we compare the results of the two 
numerical methods and the theoretical predictions of Sec. |T|. 



A. Numerical Integration Results 



To solve Eq. ([TTD numerically, we represent the field as an L-vector in which each 
element is allowed to vary continuously. For the numerical solutions presented in this work 
we have chosen L=100. (Doubling L in this case produced no discernible change in the 
results for the fields and temperatures we studied.) We set the interaction constants J=l/2 
and J'=0, so that /3c=l, and we obtained numerical results for temperatures in the range 
0.1Tc<T<0.8Tc. Setting H>0, we determine the magnetization density (pcore at the core of 
the droplet by solving /cw(0core)=/cw(0ms) under the constraint 0core>0ms- We then force 
the type-II solution by setting 4>i=(f)core, and 4>L/2=4>ras as boundary conditions. The layers 
between L/2 and L are given by the symmetry imposed by periodic boundary conditions. 
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+ . (55) 



The methods chosen to solve Eq. ([TTD are the well-known shooting and relaxation meth- 
ods. (See, e.g. Ref. [^.) Taking the explicit form of Eq. (P) for few, "we use Eq. ( |TlD 
directly to shoot from i=l, 

= - H/J + {(3J)-^ tanh-1 (/.^ , (54) 

stopping before the solution becomes catastrophic. The remainder of the system is given 
the value 0nis, and a relaxation method is employed, using a variant of Eq. ([5^), 

0j = tanh 

The values are then relaxed iteratively in a checkerboard fashion, each layer being adjusted 
by an increment A0j, until X^i |A0i|<lO~^ '^'^^ iteration. The free-energy cost AF of 
forming the droplet is then obtained using Eq. (||): 

AF ^ r -1 

^ = [i(<^i+l + /cw(0i) - /cw(0ms)J • (56) 

In Fig. H the numerical solutions for various temperatures are shown as functions of 
\\\/ Hs = {Hs — H)/Hs and are compared with the theoretical result for AF from Eq. (PU]). 
Corrections to the 0^ field theory are seen to be substantial, especially for low temperatures, 
except for a region very close to the spinodal. We can account for this rapid departure of 
the numerical solutions from the (j)^ field theory by considering the leading field- dependent 
corrections to the cubic potential used in Eq. (|T9|) . We can write the first correction to the 
integrand of Eq. (O) as jev"^, where 



'=6 



1 d^fcw 




(57) 



To leading order, the correction to the LHS of Eq. (pT]) is then ev^, and Eq. ( 126|) may be 
rewritten 

AJ^ = -A0S - ftt^^o - fe^o - + ie(2i;oM3 + m^)J dr . (58) 

The length scale of the droplet is still given approximately by ^r, so performing the inte- 
gration in the last term introduces a factor of ^r- The first correction to AF thus has the 
temperature and field dependence of ewQ^r~e|tt|~^/^| A|'^/'*. 

We can verify the presence and form of this correction by studying the behavior of 
the numerically obtained AF as |A|— >0. For a given f3, we generate a sequence of field 
deviations |An|=Ao/'n.^ (not to be confused with the transfer-matrix eigenvalues), where Aq 
is a constant, for an increasing sequence of integers n. We then calculate the free-energy 
cost of the critical droplet AF„, using the numerical integration procedure outhned above, 
for a field H = H^ — A„. If we define two sequences Pn=A.Fn/AF, where AF is the analytic 
result of Eq. (|29|) at A=A„ with no corrections assumed, and Q„ = n{Pn — 1), then the 
dominant contributions to AF„ should give the following expansions: 

P„ = 1 + Ai{p)n-^ + A2{P)n-^ + ■■■ 

Qn = A,{P) + A2{f3)n-' + ■ ■ ■ , (59) 
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where the amphtudes Ai{P) are functions only of P, and in particular, Ai(P)oce\a\~^/'^. 
For even n, we transform the numerical sequences P„ and Qn in order to remove the first 
corrections. The method we use, a variant of Neville-Aitken extrapolation |^ , is described 



in detail in Ref. |Q, where it was used to extrapolate with much success finite-size data 
carrying very substantial corrections. We apply it here to remove 0{n^^) corrections to Pn 
and Qn- 

Pn ~ 2-Pn. — Pn/2 

= 2Q„ - Q„/2 , (60) 

which give = 1 + 0(n^^) and Q'^ = Ai{P) + 0{n~'^). The constancy of the leading terms 
in Pn and Qn was verified by calculating estimates and to the exponents of n in 
these terms, which we expect to be zero: 

iy„^ = (ln2)-Mn(P;/P^/2), (61) 

and similarly for Qn, which gives the leading exponent to 0{n^'^). For each temperature, 
we computed the sequences AP„ using Ao=0.5ifs and ?t,=1, 2, 3, . . .. The sequences were 
terminated when P„ came within 0.5% of unity. The resulting estimates of the |A|-exponent 
(^rf+f) leading term lay between 1.249 and 1.251, and the estimates for the |A|- 

exponent in the first correction term (|W^^+|) lay between 1.74 and 1.75, both confirming 
the field-theoretical predictions of 5/4 and 7/4, respectively ||5^. We used the value of P^ as 
an estimate for the leading coefficient to Pn- These estimates lay between 0.9998 and 1.0001, 
confirming Eq. (^) to high accuracy. Dividing the predicted temperature dependence of 
Ai from Q'j^, we found the resulting coefficient to vary only about 15% over the entire range 
of temperatures, giving evidence that the predicted temperature dependence of the first 
correction to AF is also as predicted. For H/Hg-^l, the logarithmic cusp singularity in AF 
predicted in Refs. was also observed in our numerical integration results |5^. Since 



transfer-matrix results were not obtainable at sufficiently low fields, for this work we will 
not discuss this field region further. 



B. Constrained- Transfer-Matrix Results 

The CTM method was applied to QIDI systems of finite cross section and infinite 
length, with interaction constants J=l/2 and J'=0, so that /5c=l. The temperatures studied 
using the CTM method were in the range 0.2Tc<T<0.8Tc, and data at each temperature 
was collected for fields ranging from H=0 to H^1.2Hs. The matrix was tridiagonalized 
by a Householder reduction, and then diagonalized by a QL algorithm with implicit shifts 
||55|| . The computations were performed in 128-bit precision on a Cray Y-MP 4/32 using 
about 200 CPU-hours. Computations were limited mainly by machine precision, since the 
process described in Sec. |^ involves the reweighting and resumming of eigenvalues whose 
magnitudes range over thirty or more decades. The lower limit on the temperatures studied 
was due to limited machine precision, whereas the upper limit was due to the extremely 
large cross sections required to minimize finite-size effects. For most temperatures, systems 
were studied for up to the range 100<A^<200, although for the highest temperatures we 
were able to obtain results for A^^500. 
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In Fig. 1^ typical eigenvalue and layer magnetization spectra of the transfer matrix for 
the QIDI model are shown as functions of H . The eigenvalue spectrum is shown on a scale 
so that its lowest branch is the equilibrium free-energy density of the system. This spec- 
trum can be divided into three regions p9|. The first is a region of low- lying states with 
positive magnetization. Above this region is a region of alternately polarized states that 
become nearly degenerate at regular intervals in H. Since all the eigenvectors of the re- 
duced transfer matrix belong to the same symmetry class, the eigenvalue branches are never 
actually degenerate and therefore do not cross, but rather quickly reverse their direction. 
In Ref. the gaps between these nearly degenerate eigenvalues were found to approach 
zero exponentially with A^. The quick reversal of the eigenvalue branches can be seen in 
the near discontinuities in the magnetization branches, in which the values of the magneti- 
zations associated with the two nearly degenerate eigenvalues are effectively traded. This 
region terminates approximately at the Curie- Weiss spinodal field H^. of Eq. (|T6|). The third 
region, above the second, consists of eigenvalues that have little polarization. For the re- 
sults presented here we are concerned with the lowest-lying states of the second region, that 
is, the lowest-lying states with magnetization opposite the applied magnetic field. These 
states form a composite branch that agrees well with the values shown for the free-energy 
density and magnetization of the metastable state for the Curie- Weiss ferromagnet in the 
thermodynamic limit. This agreement improves as N is increased. 

The real and imaginary parts of fa from the CTM, as computed by Eq. are shown in 
Fig. ^ for a single system size at fixed temperature. The real part of fa exhibits a composite 
branch that is nearly identical to the composite "metastable" branch shown in Fig. ^, and 
thus is also in good agreement with the known free-energy density of the metastable state. 
In addition, the corresponding values of |Im/Q,| for this branch are nearly zero except in the 
region of the spinodal. Figure ^ shows the same branches of llm/d as Fig. ^ on a logarithmic 
scale. Note that the values of |Im/Q,| range over more than 30 decades. The structure of the 
"metastable" branch as a succession of lobes rather than as a smooth function is due to the 
mixing of nearly-degenerate eigenvectors (see Fig. Since we were interested in obtaining 
values from branches representing the single "metastable" phase we selected values of 
H for which Aq,=(Aq,+iAq,_i)^/^, thus ensuring a "safe" distance from the near-degeneracies. 
The resulting points are shown on the figure and together roughly give tangent points for a 
reasonable envelope for the lobes. We identify this envelope with Im/ms, the imaginary part 
of the constrained free-energy density for the metastable state. The envelope for a larger 
system at the same temperature is also shown. By inspecting the envelopes for many system 
sizes, we observed roughly an exponential dependence on well inside the spinodal and a 
crossover to a slower scaling near the spinodal. The evidence for this scaling behavior will 
be made clear and quantitative in the next subsection. 

The presence of the crossover can be seen clearly in Fig. where the second derivative 
of |Im/ms| with respect to H is plotted versus reduced field H/Hg for a fixed temperature 
and various values of A^. Each derivative was calculated by the midpoint method, using 
only the "tangent" points chosen as described above. As A^ increases, a singularity in 
this derivative develops and pushes closer to the spinodal. Although the value of Hs is 
temperature dependent, we observed the same behavior in |Im/ins| for the entire range of 
temperatures studied. 
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C. Finite-Range-Scaling Results 



First we turn our attention to the scaling behavior of |Im/jns| for fields in the range 
0<H<Hs, and we consider fields far enough from the spinodal that the Boltzmann weight 
g-/3AF gg^g ^YiQ scale for the metastable lifetime. Since the transfer-matrix data lack the 
extremely high resolution required to determine the prefactor to this weight, we concentrate 
on the quantity In |Im/ms|. Assuming a form for |Im/ms| as given by Eqs. (^Qf ) and (]^), we 
write 

In |Im/^s,iv(T, H)\ ~ -f3N''A{T, H) , (62) 

where the exponent a and the function A{T,H) are undetermined. 

For each system size N with A^=0(mod 4), we calculate finite-range estimates for a at 
fixed T and H assuming that the dominant correction is from the prefactor and is therefore 
0{\nNy. 

n o\~U / ln|Im/ins,jv| -ln|Im/„,s,;v/2| \ .„„x 
cT7v= ln2 In ■ — — ^ , (63) 



which gives aj\f = a + 0{N ^) [^. Whereas the values for In |Im/ins,Ar| in Eq. (|63|) were 



taken directly from the tangent points, the values of In |Im/jns| for smaller systems had to be 
interpolated to the same field value since both the number and positions of the lobes are A^- 
dependent. We performed the interpolation by fitting a quadratic form exactly to the three 
nearest points. Figure |^ shows utv for the largest values of N numerically attainable, plotted 
versus reduced field for various temperatures. Clearly these estimates are quite consistent 
with the A^-dependence (cr=l) of AF given by Eq. (pOf). Note that as the estimates approach 
Hs, they drop significantly, indicating the rapidly decreasing importance of the A^-dependent 
exponential factor of |Im/ins,Ar| compared to the algebraic prefactor, as was pointed out in 
Sec. m. 

For even A^, finite- A^ estimates for A{T,H) were calculated assuming that o"=l, and 
again that the dominant correction is O (In A^) : 

= -;^(ln |Im/ms,Ar| - In |Im/nis,7V/2|) • (64) 

Values for In |Im/ms tvI were taken as described in the previous paragraph. Figure || shows 
a typical set of estimates for AF. Due to logarithmic corrections, the convergence of naive 
estimates using Eq. (|6^) directly, A(T, if)^(/5A^)~^ In |Im/ms|, is quite slow. Using Eq. 
(|6^) gives greatly improved estimates, which agree much better with the free-energy cost 
of nucleation obtained by numerical integration, and the range of fields over which there 
is good agreement extends much closer to the spinodal. In addition, these extrapolated 
estimates show a trend near the spinodal toward the numerical interation result, although 
the limited accuracy in the interpolations used precluded any further extrapolation. 

Figure ^ shows Atv for the largest values of A^ numerically attainable, plotted versus 
reduced field for various temperatures. These are compared with the free-energy cost of 
nucleation obtained though Eq. ( |5B| ) by numerically integrating the Euler-Lagrange equation 
as described in Sec. [V A| . As can be seen from Figs. ^ and the agreement between the 



extrapolated CTM estimates and the exact results is impressive. The relative error for 
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these estimates were consistently below 2% over most of the field and temperature range 
and remained below 5% until the crossover region was entered. By comparing Fig. I with 
Fig. 0, where the free-energy cost of nucleation obtained by numerical integration is shown 
together with the 0^ field-theoretical results, one sees that the CTM estimates for AF 
faithfully reproduce the higher-order corrections to the 0^ field theory over a wide range of 
fields and temperatures. This strongly indicates the consistency between the CTM method 
and the droplet theory of nucleation, at least whenever a substantial free-energy barrier 
against nucleation exists. 

Next we consider the behavior of |Im/nis| near the spinodal, where the extrapolated 
CTM estimates A^r lie significantly above the infinite- iV free-energy cost AF obtained by 



numerical integration. As discussed in Sec. |T|, here we must carefully consider the behavior 



of the prefactor to the Boltzmann weight. Since the mean-field spinodal is a line of critical 
points in a (p^ field theory, we can apply critical finite-range scaling ||3^ in this region. By 
recasting Eqs. (|27|) , (129|) , and (|40|) in terms of a scaling variable (=R'^^\X\, we determine 
the finite-range scaling relation for |Im/| to be 

|Im/|=i?^'^'$(AC), (65) 
where the scaling function $ depends on the sign of X = H — Hg-. 



^^'^^ Ui«r^/'c^/' (i^i>i^si) 



(66) 



Note that as C^O from either direction, $ goes to zero, with a cusp singularity at |if | = |iJs| • 
This singular scaling function cannot be reproduced for finite N. We therefore expect that 
the finite- iV CTM estimates for |Im/ms| near Hs should be dominated by corrections to 
scaling. Indeed, the transfer- matrix results for |Im/ms| only showed qualitative agreement 
with the field theory for Hk-H^. Power-law scaling was found at ifg, but the exponent 
of Rn varied between roughly —0.8 at low temperatures and —1 at high temperatures, in 
significant disagreement with Eq. (|65|), and large fiuctuations in the data precluded the use 
of extrapolation techniques. In addition, the values of the scaling functions calculated from 
|Im/ins| were consistently greater than those of the theoretical scaling functions, varying 
from a factor of roughly 2 at low temperatures to 5 at high temperatures. Several possible 
explanations for these results will be considered in the next section. 



VI. DISCUSSION 

In this work we have applied a scalar field theory and the recently developed constrained- 
transfer-matrix (CTM) method to the study of metastability in the ferromagnetic quasi-one- 
dimensional Ising model, a two-state model with weak, long-range interactions. By solving 
the Euler-Lagrange equation associated with the Ginzburg-Landau- Wilson Hamiltonian for 
the model, we have identified the equilibrium and metastable states as spatially uniform 
configurations with the order parameter taking the Curie- Weiss mean-field values, and we 
have identified the critical fiuctuation through which the metastable state decays. By com- 
puting the free-energy cost of the critical fiuctuation numerically, we were able to calculate 
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the "Boltzmann weight" that appears in the analytic continuation of the free energy across 
the first-order phase transition for the entire region of metastabihty. Except in the region of 
the classical spinodal, this weight gives the dominant time scale for decay of the metastable 
state. In the region of the spinodal we have mapped the Hamiltonian to a 0'^ field theory 
to obtain an expression for the analytically continued free energy with no undetermined 
parameters. Assuming the applicability of the Fokker-Planck equation to the dynamics, we 
have thus also obtained an exact expression for the decay rate of the metastable state for 
this model near the classical spinodal. In performing the analytic continuation of the free 
energy, we have found that its corrections are quite substantial, leading to large differences 
in, for example, the free-energy cost of nucleation as one moves away from the spinodal to- 
ward the first-order transition. Numerical solution of the Euler-Lagrange equation strongly 
supports the field-theoretical result for the free energy-cost of nucleation, including the first 
correction term. 

We have outlined a method by which a complex analogue of the free energy for a con- 
strained system is obtained directly from the transfer matrix for the unconstrained system. 
Extensive symmetry reductions have enabled us to collect transfer-matrix data for systems 
with very large cross sections, and hence very large interaction ranges. We found that 
the real part of the constrained free-energy density associated with the metastable phase, 
Re/ms, rapidly approaches the free-energy density of the metastable phase as N—>-oo. The 
associated imaginary part |Im/jns| is extremely small, showing exponential dependence on 
the interaction range over most of the region of metastabihty and a crossover to power-law 
scaling near the classical spinodal. We found strong evidence for this crossover in the be- 
havior of the second field derivative of |Im/ms|- Using numerical extrapolation techniques, 
we have demonstrated that over a wide range of fields and temperatures within the region 
of metastabihty, the complex "constrained" free-energy density obtained by this method 
agrees very well with the behavior of the analytically continued free-energy density. The 
estimated free-energy cost of the critical droplet was found to lie within 2% of the value 
obtained by solving the Euler-Lagrange equation numerically. In fact, for low temperatures 
and for fields sufficiently far from the spinodal, the transfer-matrix data were found to agree 
with the numerical solution much more closely than the analytic field-theoretical results. 

However, we found that the constrained free-energy density, as defined by Eq. using 
the reweighting scheme defined in Eq. (^Tj), does not show consistent finite-range scaling at 
the spinodal. There are several possible explanations for this result. If we consider the 
possibility that finite-size corrections to the free-energy cost AF of nucleation are affecting 
the exponents in the prefactor of |Im/|, then those corrections must be O(lnA^) |5^. (This 
rules out, for example, the effect of relaxing Stirling's approximation, which was used to 
derive Eq. (^, since it gives only a correction of order unity.) Another possibility is that 
the effects of eigenvector mixing are too strong in this region to obtain reasonable estimates 
for |Im/nis|. In addition, for finite cross sections N, the value of the spinodal may not be 
entirely real, so that a dominant contribution to llm/msl emerges from the CTM data due 
to finite-size rounding. If f^s does represent the analytic continuation of the free-energy 
density, then it is possible that the eigenvalues are being reweighted incorrectly in this 
region, allowing unwanted or unphysical fiuctuations. These and other possibilities will be 
explored in further work. 

The CTM method clearly shows promise in the characterization of metastable phases. As 
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a nonperturbative method, it treats all possible fluctuations in a single calculation, and hence 
does not require any particular assumptions about the detailed structure of the partition 
function. Since it automatically identifles the fluctuations that are important to nucleation, 
the CTM method has potential applications to the study of metastability in, for example, 
disordered systems, in which critical fluctuations are difficult to characterize. The usefulness 
of this method has been demonstrated in a study of the 2D Ising ferromagnet ||3^ , in which 
excellent quantitative agreement was found between fleld-theoretical and CTM estimates of 
the surface free energy and the shape of the critical droplets, as well as consistency with 
Monte-Carlo estimates of the metastable lifetime Furthermore, in a three-state 

model with weak, long-range interactions |Q, |Im/nis| appears to remain consistent with 
the metastable lifetime, even under conditions of competing metastable states, where it has 
been argued that the analytic continuation of the free energy may no longer be a valid 
measure of the lifetime. A Monte-Carlo study of metastability in the QIDI model, using 
recently developed techniques, is in progress, and a further investigation of the connection 
between the dynamics of metastable states and the transfer matrix is planned. 



ACKNOWLEDGMENTS 



We would like to thank T. Fiig and C. C. A. Giinther for useful discussions. This work 
was supported in part by Florida State University through the Supercomputer Computations 
Research Institute (Department of Energy Contract No. DE-FC05-85ER25000), through the 
Center for Materials Research and Technology, and through Cray Y-MP supercomputer time, 
and by National Science Foundation Grant No. DMR-9013107. 



21 



REFERENCES 



[1] J. D. Gunton, M. San Miguel, and P. S. Sahni, in Phase Transitions and Critical Phe- 
nomena, Vol. 8, edited by C. Domb and J. L. Lebowitz (Academic, New York, 1983). 
[2] K. Binder, Rep. Prog. Phys. 50, 783 (1987). 
[3] J. S. Langer, Ann. Phys. 41, 108 (1967). 
[4] J. S. Langer, Phys. Rev. Lett. 21, 973 (1968). 
[5] J. S. Langer, Ann. Phys. 54, 258 (1969). 
[6] K. Binder and E. StoU, Phys. Rev. Lett. 31, 47 (1973). 
[7] K. Binder, Phys. Rev. B 8, 3423 (1973). 

[8] K. Binder and H. Miiller-Krumbhaar, Phys. Rev. B 9, 2328 (1974). 
[9] R. J. McCraw and L. S. Schulman, J. Stat. Phys. 18, 293 (1978). 
[10] R. J. McCraw, Phys. Lett. 75A, 379 (1980). 

[11] C. M. Newman and L. S. Schulman, J. Stat. Phys. 23, 131 (1980). 

[12] G. Rocpstorff and L. S. Schulman, J. Stat. Phys. 34, 35 (1984). 

[13] M. Biittiker and R. Landauer, Phys. Rev. Lett. 43, 1453 (1979). 

[14] M. Biittiker and R. Landauer, Phys. Rev. A 23, 1397 (1981). 

[15] W. Klein and G. Unger, Phys. Rev. B 28, 445 (1983). 

[16] G. Unger and W. Klein, Phys. Rev. B 29, 2698 (1984). 

[17] G. G. Gallan and S. Goleman, Phys. Rev. D 16, 1762 (1977). 

[18] N. J. Giinther, D. A. Nicole, and D. J. Wallace, J. Phys. A 13, 1755 (1980). 

[19] W. N. Gottingham, D. Kalafatis, and R. Vinh Mau, Phys. Rev. B 48, 6788 (1993). 

[20] R. J. Baxter and I. G. Enting, J. Stat. Phys. 21, 103 (1979). 

[21] G. A. Baker, Jr. and D. Kim, J. Phys. A 13, L103 (1980). 

[22] M. J. Lowe and D. J. Wallace, J. Phys. A 13, L381 (1980). 

[23] 1. G. Enting and R. J. Baxter, J. Phys. A 13, 3723 (1980). 

[24] G. J. Hamer, J. Phys. A 16, 1247 (1983). 

[25] B. Gaveau and L. S. Schulman, Lett. Math. Phys. 18, 201 (1989). 

[26] G. M. Newman and L. S. Schulman, J. Math. Phys. 18, 23 (1977). 

[27] V. Privman and L. S. Schulman, J. Phys. A 15, L231 (1982). 

[28] V. Privman and L. S. Schulman, J. Stat. Phys. 31, 205 (1982). 

[29] M. A. Novotny W. Klein, and P. A. Rikvold, Phys. Rev. B 33, 7729 (1986). 

[30] P. A. Rikvold, Prog. Theor. Phys. Suppl. 99, 95 (1989). 

[31] P. A. Rikvold, Phys. Scr. T38, 36 (1991). 

[32] P. A. Rikvold, B. M. Gorman, and M. A. Novotny AIP Gonf. Proc. 256, 549 (1992). 
[33] G. G. A. Giinther, P. A. Rikvold, and M. A. Novotny, Phys. Rev. Lett., 13 Dec. (1993). 
[34] P. A. Rikvold, B. M. Gorman, and M. A. Novotny Phys. Rev. E 47, 1474 (1993). 
[35] H. Wagner and H. Horner, Adv. Phys. 23, 587 (1974). 
[36] A. J. Leggett, Rev. Mod. Phys. 47, 331 (1975). 

[37] P. G. de Gennes, Scaling Concepts in Polymer Physics (Gornell Univ. Press, Ithaca, 
1979). 

[38] K. Binder, Phys. Rev. A 29, 341 (1984). 
[39] R. A. Gowley, Adv. Phys. 29, 1 (1980). 

[40] J. R. Schrieffer, Theory of Superconductivity (Addison- Wesley, New York, 1983). 
[41] P. B. Littlcwood and P. Ghandra, Phys. Rev. Lett. 57, 2415 (1986). 
[42] P. Ghandra, Phys. Rev. A 39, 3672 (1989). 



22 



R. B. Griffiths, C.-Y. Wang, and J. S. Langer, Phys. Rev. 149, 301 (1966). 
W. Paul, D. W. Heermann, and K. Binder, J. Phys. A 22, 3325 (1989). 
J. D. Gunton and M. C. Yalabik, Phys. Rev. B 18, 6199 (1978). 
T. S. Ray and W. Klein, J. Stat. Phys. 61, 891 (1990). 
T. S. Ray, J. Stat. Phys. 62, 463 (1991). 

L. S. Schulman, in Techniques and Applications of Path Integration (Wiley, New York, 
1981), Chap. 29. 

The expression for vi was misprinted in Eq. A2 of Ref. ||14|| . 

The results of Ref. [|TB| should be corrected as follows: (1) the exponent 1/4 in Eq. 4.8 
should be —1/4; (2) the exponent 3(i/8 in Eq. 4.11 should be d{Q—d)/9>] and consequently 
(3) the exponent (l/2+5rf/8) in Eq. 4.16 should be rf(l-rf/8). 
C. Domb, Adv. Phys. 9, 149 (1960). 

M. Marcus and H. Mine, in A Survey of Matrix Theory and Matrix Inequalities (Dover, 
New York, 1964), p. 124. 

R. E. Blahut, in Principles and Practice of Information Theory (Addison- Wesley, Read- 
ing, 1987), pp. 61-64. 

T. Fiig, B. M. Gorman, P. A. Rikvold, and M. A. Novotny (FSU-SCRI preprint). 

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes 

(Cambridge Univ. Press, New York, 1992). 

A. J. Guttmann, in Phase Transitions and Critical Phenomena, Vol. 13, edited by C. 
Domb and J. L. Lebowitz (Academic, New York, 1990). 

B. M. Gorman (unpublished). 

P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides (FSU-SCRI preprint). 

C. C. A. Giinther, P. A. Rikvold, and M. A. Novotny (unpubhshed) . 



23 



FIGURES 



FIG. 1. Schematic diagrams of the solutions to Eq. (14) for the cases of H=0 and H>0. The 
solutions are shown together with sketches of the corresponding Curie- Weiss free-energy density 
/cw- The diamonds mark the uniform (type-I) solutions. The dashed horizontal line in each case 
represents the interval over which the order parameter ranges in the type-II solution. The hashed 
regions represent the bands of type-Ill solutions, which oscillate in space between magnetization 
densities with the same value of /cw- Below these sketches of the intervals are sketches in real 
space of the type-II solutions. 



FIG. 2. Estimates for the free-energy cost AF of forming the critical droplet in the QIDI 
model, shown as a cross-section {N-) density and as a function of the H-field deviation from the 
spinodal, \X\/Hs = (Hg — H)/Hs. The solid curves are the results of numerically solving Eq. (11) 
at temperatures T/Tc=0.2, 0.4, 0.6, and 0.8. The curves corresponding to this order run from top 
to bottom. The dashed straight lines show the estimates from the 0^ field-theoretical result of Eq. 
( p9| ) at the same temperatures. These lines appear in the same order, with the two top lines nearly 
coincident. 



FIG. 3. Eigenvalue (top) and magnetization (bottom) spectra computed from the transfer 
matrix for a 35xoo QIDI system at T=0.5Tc and plotted versus H. The value of the Curie- Weiss 
spinodal field Hg is marked by the dotted vertical line, and the values of the free-energy density and 
magnetization for the equilibrium, metastable, and unstable stationary states of the Curie- Weiss 
ferromagnet are marked by the thick gray lines. Fifteen eigenvalue branches lie above the chosen 
plotting range, so they are not visible. See Sec. VB for a detailed description. 



FIG. 4. All branches of Rc/q, (top) and llm/d (bottom) computed from constrained transfer 
matrices for a 35xoo QIDI system at T=0.5Tc and plotted versus H. The dotted vertical line 
indicates the spinodal field. The branches of Re/a are plotted on the same scale as the eigenvalue 
branches of Fig. ^, so that the two spectra may be compared directly. 



FIG. 5. Branches of ln|Im/a| computed from constrained transfer matrices for a 35xoo QIDI 
system at T=0.5Tc and plotted versus H. The (-I-) symbols touching the lobes are calculated 
at fields chosen by the criterion AQ,=(Ao,+iAa_i)^/^ for the "metastable" eigenvector \a). These 
points are joined by straight lines. A second set of points, for N=70, is shown to illustrate the 
A^-scaling behavior of ln|Im/Q,|. This second set is expected to continue downward at lower fields 
roughly along the dashed line, but the numerical precision was not sufficient to resolve points in 
the low-field region. The spinodal field is marked by the dotted vertical line. 



FIG. 6. The second derivative of llm/d with respect to H for QIDI systems with various cross 
sections N at r=0.2rc. The symbols are the results of the midpoint method applied to the points 
used to define the envelope over the lobes of |Im/Q,|. (See Fig. |5[) The lines joining these points 
are splines serving only as guides to the eye. The dotted vertical line indicates the spinodal field. 
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FIG. 7. Finite- estimates for the exponent a in Eq. 
against reduced field H/Hs for various temperatures T. 
field-theoretical value (o"=l) over this range. 



(|62|), for the largest A'' attainable, plotted 
The dotted horizontal line indicates the 



FIG. 8. Finite-A^ estimates for A{T,H) in Eq. ( |62D with a=l fixed, plotted against 
|A|/ifs = {Hs — H)/Hs for T=0.5Tc. Each thin line connects the naive estimates one would obtain 
for a particular N directly from Eq. (|62[). These estimates decrease in value as N increases. The 
data points are estimates using Eq. (|6j). The diamonds are results using the largest attainable. 
The circles are results for half of the largest A^ and are plotted only in the region where finite-size 
effects are large. The thick line indicates the free-energy cost of forming the numerical "critical 
droplet" solution of Eq. (Ill]), obtained by numerical integration as described in Sec. V A. 



FIG. 9. Finite-A^ estimates for A{T,H) in Eq. (H) with a=l fixed, extrapolated by Eq. (g) 
for the largest A^ attainable, and plotted against |A|/i7s = {Hg — H)/Hs for various temperatures 
T. The lines indicate the free-energy cost of forming the numerical "critical droplet" solution of 
Eq. (In]) at the same temperatures, obtained by numerical integration as described in Sec. |V A| . 
Compare with Fig. 0. 
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